# "Media's Influence on LGBTQ Support Across Africa" by Stephen Winkler

# Code to run robustness models and create:
#   - Table A.7 (logit results for religion placebo)
#   - Table A.8 (logit results for ethnicity placebo)
#   - Table A.9 (logit results for HIV+ placebo)
#   - Table A.10 (logit results for immigrant placebo)
#   - Table A.11 (ols results for freedom house interaction)
#   - Table A.12 (logit results for freedom house interaction)
#   - Table A.15 (ols results for KOF interaction)
#   - Table A.16 (logit results for KOF interaction)

rm(list=ls()) # clear environment
source("Rprofile.R") # setup Rprofile (load packages, etc.)

# Load data
myData <- readRDS("data/afrobarometer.rds")

# Explanatory variables:
myData$radio <- as.numeric(myData$radio)
myData$tv <- as.numeric(myData$tv)
myData$newspaper <- as.numeric(myData$newspaper)
myData$internet <- as.numeric(myData$internet)
myData$social_media <- as.numeric(myData$social_media)
myData$media_aggregate <- as.numeric(myData$media_aggregate)
myData$media_noradio <- as.numeric(myData$media_noradio)
myData$media_notv <- as.numeric(myData$media_notv)
myData$media_nopaper <- as.numeric(myData$media_nopaper)
myData$media_nointernet <- as.numeric(myData$media_nointernet)
myData$media_nosocialmedia <- as.numeric(myData$media_nosocialmedia)

###############################################
# Robust Check: Effect on other Tolerance DVs
# (Placebo Tests)
###############################################

# Model 1: Media Aggregate (Controls + CFE + DCSE)
# Religion
model <-  #linear model
  relig_tol2 ~
  media_aggregate +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.all.relig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.all.relig

# Ethnicity
model <- 
  ethnic_tol2 ~
  media_aggregate +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.all.ethnic <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.all.ethnic

# Homosexuality
model <- 
  sexuality2 ~
  media_aggregate +
  ethnic_tol2 + relig_tol2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.all.sexuality <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.all.sexuality

# HIV
model <- 
  hiv_tol2 ~
  media_aggregate +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.all.hiv <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.all.hiv

# Immigrants
model <- 
  immig_tol2 ~
  media_aggregate +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.all.immig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.all.immig

# Model 2: Media Types (Control other media + controls + cfe + dcse)

# Radio + Relig
model <-
  relig_tol2 ~
  radio + media_noradio +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.radio.relig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.radio.relig

# Radio + Ethnic
model <-
  ethnic_tol2 ~
  radio + media_noradio +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.radio.ethnic <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.radio.ethnic

# Radio + HIV
model <-
  hiv_tol2 ~
  radio + media_noradio +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.radio.hiv <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.radio.hiv

# Radio + Sexuality
model <-
  sexuality2 ~
  radio + media_noradio +
  ethnic_tol2 + relig_tol2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.radio.sexuality <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.radio.sexuality

# Radio + Immig
model <-
  immig_tol2 ~
  radio + media_noradio +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.radio.immig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.radio.immig

# TV + Relig
model <-
  relig_tol2 ~
  tv + media_notv +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.tv.relig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.tv.relig

# TV + Ethnic
model <-
  ethnic_tol2 ~
  tv + media_notv +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.tv.ethnic <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.tv.ethnic

# TV + HIV
model <-
  hiv_tol2 ~
  tv + media_notv +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.tv.hiv <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.tv.hiv

# TV + Sexuality
model <-
  sexuality2 ~
  tv + media_notv +
  ethnic_tol2 + relig_tol2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.tv.sexuality <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.tv.sexuality

# TV + Immig
model <-
  immig_tol2 ~
  tv + media_notv +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.tv.immig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.tv.immig

# Newspaper + Relig
model <-
  relig_tol2 ~
  newspaper + media_nopaper +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.newspaper.relig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.newspaper.relig

# Newspaper + Ethnic
model <-
  ethnic_tol2 ~
  newspaper + media_nopaper +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.newspaper.ethnic <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.newspaper.ethnic

# Newspaper + HIV
model <-
  hiv_tol2 ~
  newspaper + media_nopaper +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.newspaper.hiv <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.newspaper.hiv

# Newspaper + Sexuality
model <-
  sexuality2 ~
  newspaper + media_nopaper +
  ethnic_tol2 + relig_tol2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.newspaper.sexuality <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.newspaper.sexuality

# Newspaper + Immig
model <-
  immig_tol2 ~
  newspaper + media_nopaper +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.newspaper.immig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.newspaper.immig

# Internet + Relig
model <-
  relig_tol2 ~
  internet + media_nointernet +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.internet.relig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.internet.relig

# Internet + Ethnic
model <-
  ethnic_tol2 ~
  internet + media_nointernet +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.internet.ethnic <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.internet.ethnic

# Internet + HIV
model <-
  hiv_tol2 ~
  internet + media_nointernet +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.internet.hiv <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.internet.hiv

# Internet + Sexuality
model <-
  sexuality2 ~
  internet + media_nointernet +
  ethnic_tol2 + relig_tol2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.internet.sexuality <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.internet.sexuality

# Internet + Immig
model <-
  immig_tol2 ~
  internet + media_nointernet +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.internet.immig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.internet.immig

# Soc Media + Relig
model <-
  relig_tol2 ~
  social_media + media_nosocialmedia +
  ethnic_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.social_media.relig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.social_media.relig

# Soc Media + Ethnic
model <-
  ethnic_tol2 ~
  social_media + media_nosocialmedia +
  relig_tol2 + sexuality2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.social_media.ethnic <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.social_media.ethnic

# Soc Media + HIV
model <-
  hiv_tol2 ~
  social_media + media_nosocialmedia +
  ethnic_tol2 + sexuality2 + relig_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.social_media.hiv <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.social_media.hiv

# Soc Media + Sexuality
model <-
  sexuality2 ~
  social_media + media_nosocialmedia +
  ethnic_tol2 + relig_tol2 + hiv_tol2 + immig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.social_media.sexuality <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.social_media.sexuality

# Soc Media + Immig
model <-
  immig_tol2 ~
  social_media + media_nosocialmedia +
  ethnic_tol2 + sexuality2 + hiv_tol2 + relig_tol2 + # social tolerance of other groups
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
lm.result <- lm(model, data = mdata)
lm.social_media.immig <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.social_media.immig


# Generate Latex Tables
# Religion as DV: Table A.7
stargazer(lm.all.relig, lm.radio.relig, lm.tv.relig, lm.newspaper.relig, lm.internet.relig, lm.social_media.relig,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          title = "Effect of Media Consumption on Religious Tolerance (OLS)",
          omit= "country",
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level."))
# Ethnicity as DV: Table A.8
stargazer(lm.all.ethnic, lm.radio.ethnic, lm.tv.ethnic, lm.newspaper.ethnic, lm.internet.ethnic, lm.social_media.ethnic,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          title = "Effect of Media Consumption on Ethnic Tolerance (OLS)",
          omit= "country",
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level."))
# HIV as DV: Table A.9
stargazer(lm.all.hiv, lm.radio.hiv, lm.tv.hiv, lm.newspaper.hiv, lm.internet.hiv, lm.social_media.hiv,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          title = "Effect of Media Consumption on HIV+ Tolerance (OLS)",
          omit= "country",
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level."))
# Immig as DV: Table A.10
stargazer(lm.all.immig, lm.radio.immig, lm.tv.immig, lm.newspaper.immig, lm.internet.immig, lm.social_media.immig,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          #column.labels=c("Logit", ""), # "Pol Herf"),
          title = "Effect of Media Consumption on Immigrant Tolerance (OLS)",
          omit= "country",
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level."))

############################################
# Robustness Check : Free Press Interaction
############################################
# Rerun main models but with an interaction for free press

# All Media
model <-
  sexuality2 ~
  media_aggregate + media_aggregate:fh_scale +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.all <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.all
lm.result <- lm(model, data = mdata)
lm.result.robust.all <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.all

# Radio
model <-
  sexuality2 ~
  radio + radio:fh_scale  +
  media_noradio +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.radio <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.radio
lm.result <- lm(model, data = mdata)
lm.result.robust.radio <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.radio

# TV
model <-
  sexuality2 ~
  tv + tv:fh_scale  +
  media_notv +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.tv <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.tv
lm.result <- lm(model, data = mdata)
lm.result.robust.tv <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.tv

# Newspaper
model <-
  sexuality2 ~
  newspaper + newspaper:fh_scale  +
  media_nopaper +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.newspaper <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.newspaper
lm.result <- lm(model, data = mdata)
lm.result.robust.newspaper <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.newspaper

# Internet
model <-
  sexuality2 ~
  internet + internet:fh_scale  +
  tolerance_noLGBT +
  media_nointernet +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.internet <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.internet
lm.result <- lm(model, data = mdata)
lm.result.robust.internet <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.internet

# Social Media
model <-
  sexuality2 ~
  social_media + social_media:fh_scale  +
  media_nosocialmedia +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.social_media <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.social_media
lm.result <- lm(model, data = mdata)
lm.result.robust.social_media <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.social_media

# Generate OLS Latex Table: A.11
stargazer(lm.result.robust.all, lm.result.robust.radio, lm.result.robust.tv, lm.result.robust.newspaper,
          lm.result.robust.internet, lm.result.robust.social_media,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          title = "Effect of Media Consumption on LGBT Attitudes (OLS Models)",
          omit= "country",
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level.")
)

# Generate Logit Latex Table: A.12
stargazer(logit.result.robust.all, logit.result.robust.radio, logit.result.robust.tv, logit.result.robust.newspaper,
          logit.result.robust.internet, logit.result.robust.social_media,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          #column.labels=c("Logit", ""), # "Pol Herf"),
          title = "Effect of Media Consumption on LGBT Attitudes (Logit Models)",
          omit= "country",
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level.")
)


############################################
# Robustness Check: KOF Interaction
############################################
# Rerun main models but with an interaction for globalization

# All Media
model <-
  sexuality2 ~
  media_aggregate + media_aggregate:KOFSoGI +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.all <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.all
lm.result <- lm(model, data = mdata)
lm.result.robust.all <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.all

# Radio
model <-
  sexuality2 ~
  radio + radio:KOFSoGI  +
  media_noradio +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.radio <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.radio
lm.result <- lm(model, data = mdata)
lm.result.robust.radio <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.radio

# TV
model <-
  sexuality2 ~
  tv + tv:KOFSoGI  +
  media_notv +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.tv <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.tv
lm.result <- lm(model, data = mdata)
lm.result.robust.tv <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.tv

# Newspaper
model <-
  sexuality2 ~
  newspaper + newspaper:KOFSoGI  +
  media_nopaper +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.newspaper <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.newspaper
lm.result <- lm(model, data = mdata)
lm.result.robust.newspaper <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.newspaper

# Internet
model <-
  sexuality2 ~
  internet + internet:KOFSoGI  +
  tolerance_noLGBT +
  media_nointernet +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.internet <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.internet
lm.result <- lm(model, data = mdata)
lm.result.robust.internet <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.internet

# Social Media
model <-
  sexuality2 ~
  social_media + social_media:KOFSoGI  +
  media_nosocialmedia +
  tolerance_noLGBT +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) 
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno for reference
logit.result <- glm(model, family = binomial, data = mdata)
logit.result.robust.social_media <- coeftest(logit.result, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.robust.social_media
lm.result <- lm(model, data = mdata)
lm.result.robust.social_media <- coeftest(lm.result, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.robust.social_media

# Generate OLS Latex Table: A.15
stargazer(lm.result.robust.all, lm.result.robust.radio, lm.result.robust.tv, lm.result.robust.newspaper,
          lm.result.robust.internet, lm.result.robust.social_media,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          title = "Effect of Media Consumption on LGBT Attitudes (OLS Models)",
          label = "ols_kof_table",
          omit= "country",
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level.")
)

# Generate Logit Latex Table: A.16
stargazer(logit.result.robust.all, logit.result.robust.radio, logit.result.robust.tv, logit.result.robust.newspaper,
          logit.result.robust.internet, logit.result.robust.social_media,
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          title = "Effect of Media Consumption on LGBT Attitudes (Logit Models)",
          omit= "country",
          label = "logit_kof_table",
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level.")
)